function out = RN_OptionPsi(kappa,zeta,eta,sigma,rho_sv,tau,K,r,S,v,X,alpha,phi,j)

% 
%    Returns the conditional probabilities Psi1 and Psi2

%    kappa      = speed of mean reversion stock return varaince under RN measure
%    zeta       = long-run mean coefficient 0 of stock return varaince under RN measure
%    eta       = long-run mean coefficient 1 of stock return varaince under RN measure
%    sigma      = volatility parameter of stock return variance
%    rho_sv     = correlation between stock return and its variance

%    Option Contract features
%    tau        = time to maturity 
%    K          = strike price
%    r          = riskless interest rate
%    S          = current stock price
%    v          = current stock return variance

%    X          = current sentiment
%    alpha      = weight

%    phi        = integration variable
%    j          = 1 for returning Psi1 or 2  for returning Psi2

i=sqrt(-1);
s = log(S);
if j==1
	u = 0.5;
	y = kappa - rho_sv*sigma;
else
	u = -0.5; 
	y = kappa;
end

% Defining constants in ODE for B
P=0.5*(2*u*phi*i-phi^2);
Q=y-(rho_sv*sigma*phi*i);
R=0.5*sigma^2;

P_=alpha*(u+(phi*i));
Q_=alpha*rho_sv*sigma;
R_=0.5*alpha^2;


% First-order Discretization scheme for ODEs A, B and C
h = 1/500;              % time interval choose 1/500 for the paper
t = 0:h:tau;            % time vector from 0 to maturity tau 
J=length(t);

vec_A = zeros(1,J);     % Empty vector for the deterministic function A
vec_B = zeros(1,J);     % Empty vector for the deterministic function B
vec_C = zeros(1,J);     % Empty vector for the deterministic function B
vec_A(1)=0;             % initilization
vec_B(1)=0;             % initilization
vec_C(1)=0;             % initilization

for n = 1:J-1    
    vec_C(n+1) = vec_C(n)+((eta*vec_B(n))-(alpha*vec_C(n)))*h;
    vec_B(n+1) = vec_B(n)+((P-(Q*vec_B(n))+(R*vec_B(n)*vec_B(n)))*h)...
                         +((P_*vec_C(n))+(Q_*vec_B(n)*vec_C(n))+(R_*vec_C(n)*vec_C(n)))*h;   
    vec_A(n+1) = vec_A(n)+((r*phi*i)+(zeta*vec_B(n))+(alpha*r*vec_C(n)))*h;
end
A=vec_A(end);           % Value of the deterministic function A 
B=vec_B(end);           % Value of the deterministic function B 
C=vec_C(end);           % Value of the deterministic function C

varphi = exp(A+(B*v)+(C*X)+(i*phi*s));   % characteristic function
out = real(exp(-i*phi*log(K))*varphi/(i*phi)); % returning the real part 


